SVI Mapping • findSVI

library(findSVI)
library(dplyr)
#> 
#> Attaching package: 'dplyr'
#> The following objects are masked from 'package:stats':
#> 
#>     filter, lag
#> The following objects are masked from 'package:base':
#> 
#>     intersect, setdiff, setequal, union
library(purrr)
library(leaflet)
library(knitr)
library(tidyr)
library(glue)
library(htmltools)
library(sf) #not suggests
#> Linking to GEOS 3.10.2, GDAL 3.4.1, PROJ 8.2.1; sf_use_s2() is TRUE

Get census data with geometry

First, we are using get_census_data() to obtain ZCTA-level data with simple feature geometry for PA for 2020 (Census API required).

pa_zcta_2020_geo_data <- get_census_data(year = 2020, state = "PA", geography = "zcta", geometry = TRUE)

Here, we are showing the first 10 rows of the data. With the geometry = TRUE, we’ll get a tibble with an additional column containing simple feature geometry (MULTIPOLYGON).

#> Simple feature collection with 10 features and 134 fields
#> Geometry type: MULTIPOLYGON
#> Dimension:     XY
#> Bounding box:  xmin: -80.48268 ymin: 40.10397 xmax: -79.73383 ymax: 40.83308
#> Geodetic CRS:  NAD83
#> # A tibble: 10 × 135
#>    GEOID NAME        B06009_002E B06009_002M B09001_001E B09001_001M B11012_010E
#>    <chr> <chr>             <dbl>       <dbl>       <dbl>       <dbl>       <dbl>
#>  1 15001 ZCTA5 15001        1655         417        5450         582         724
#>  2 15003 ZCTA5 15003         441         125        2352         245         250
#>  3 15004 ZCTA5 15004          62          43          63          70          21
#>  4 15005 ZCTA5 15005         293          92        1380         200          65
#>  5 15006 ZCTA5 15006          37          58          13          28           0
#>  6 15007 ZCTA5 15007           0          11          22          35           0
#>  7 15009 ZCTA5 15009         404         121        3253         384         206
#>  8 15010 ZCTA5 15010        1237         242        5636         444         719
#>  9 15012 ZCTA5 15012         638         180        2750         302         336
#> 10 15014 ZCTA5 15014         136          75         527          96         105
#> # ℹ 128 more variables: B11012_010M <dbl>, B11012_015E <dbl>,
#> #   B11012_015M <dbl>, B16005_007E <dbl>, B16005_007M <dbl>, B16005_008E <dbl>,
#> #   B16005_008M <dbl>, B16005_012E <dbl>, B16005_012M <dbl>, B16005_013E <dbl>,
#> #   B16005_013M <dbl>, B16005_017E <dbl>, B16005_017M <dbl>, B16005_018E <dbl>,
#> #   B16005_018M <dbl>, B16005_022E <dbl>, B16005_022M <dbl>, B16005_023E <dbl>,
#> #   B16005_023M <dbl>, B16005_029E <dbl>, B16005_029M <dbl>, B16005_030E <dbl>,
#> #   B16005_030M <dbl>, B16005_034E <dbl>, B16005_034M <dbl>, …

Get SVI with geometry

After getting the data ready, we can supply this tibble with simple feature geometry to get_svi().

pa_zcta_2020_geo_svi <- get_svi(2020, data = pa_zcta_2020_geo_data)

get_svi() will return the full SVI table with every SVI variables, intermediate percentile ranks, theme-specific SVIs and overall SVI (consistent with CDC/ATSDR SVI database, without MOE).

pa_zcta_2020_geo_svi %>% head(10) %>% kable()

GEOID	NAME	geometry	E_TOTPOP	E_HU	E_HH	E_POV150	E_UNEMP	E_HBURD	E_NOHSDP	E_UNINSUR	E_AGE65	E_AGE17	E_DISABL	E_SNGPNT	E_LIMENG	E_MINRTY	E_MUNIT	E_MOBILE	E_CROWD	E_NOVEH	E_GROUPQ	EP_POV150	EP_UNEMP	EP_HBURD	EP_NOHSDP	EP_UNINSUR	EP_AGE65	EP_AGE17	EP_DISABL	EP_SNGPNT	EP_LIMENG	EP_MINRTY	EP_MUNIT	EP_MOBILE	EP_CROWD	EP_NOVEH	EP_GROUPQ	EPL_POV150	EPL_UNEMP	EPL_HBURD	EPL_NOHSDP	EPL_UNINSUR	EPL_AGE65	EPL_AGE17	EPL_DISABL	EPL_SNGPNT	EPL_LIMENG	EPL_MINRTY	EPL_MUNIT	EPL_MOBILE	EPL_CROWD	EPL_NOVEH	EPL_GROUPQ	SPL_theme1	SPL_theme2	SPL_theme3	SPL_theme4	RPL_theme1	RPL_theme2	RPL_theme3	RPL_theme4	SPL_themes	RPL_themes
15001	ZCTA5 15001	MULTIPOLYGON (((-80.43758 4…	31129	16070	14093	5567	787	3037	1655	1319	6395	5450	4871	815	142	5764	955	537	94	1189	619	18.2	4.8	21.5	7.0	4.3	20.5	17.5	15.9	5.8	0.5	18.5	5.9	3.3	0.7	8.4	2.0	0.5313	0.5595	0.5174	0.4001	0.4804	0.5516	0.3258	0.5877	0.6929	0.6555	0.8173	0.7801	0.4422	0.5305	0.7071	0.8031	2.4887	2.8135	0.8173	3.2630	0.5028	0.7418	0.8173	0.8224	9.3825	0.8059
15003	ZCTA5 15003	MULTIPOLYGON (((-80.2366 40…	11212	6084	5104	2656	212	1234	441	534	2090	2352	1656	356	86	1771	196	38	8	587	26	23.7	3.5	24.2	5.5	4.8	18.6	21.0	14.8	7.0	0.8	15.8	3.2	0.6	0.2	11.5	0.2	0.7039	0.3761	0.6615	0.2798	0.5417	0.4240	0.5840	0.4991	0.7989	0.7469	0.7860	0.6724	0.2872	0.3886	0.8234	0.5585	2.5630	3.0529	0.7860	2.7301	0.5380	0.8473	0.7860	0.6515	9.1320	0.7707
15004	ZCTA5 15004	MULTIPOLYGON (((-80.38652 4…	380	207	138	75	0	6	62	38	42	63	61	21	0	30	0	25	0	0	0	19.7	0.0	4.3	22.5	10.0	11.1	16.6	16.1	15.2	0.0	7.9	0.0	12.1	0.0	0.0	0.0	0.5837	0.0000	0.0490	0.9393	0.8597	0.0959	0.2690	0.5997	0.9641	0.0000	0.6118	0.0000	0.7499	0.0000	0.0000	0.0000	2.4317	1.9287	0.6118	0.7499	0.4745	0.2520	0.6118	0.1266	5.7221	0.2287
15005	ZCTA5 15005	MULTIPOLYGON (((-80.24356 4…	9191	4275	3948	916	123	609	293	250	2534	1380	1288	92	19	459	134	23	0	135	213	10.2	2.5	15.4	4.1	2.8	27.6	15.0	14.3	2.3	0.2	5.0	3.1	0.5	0.0	3.4	2.3	0.2295	0.2593	0.2234	0.1896	0.2890	0.8678	0.1986	0.4594	0.2963	0.5312	0.4801	0.6684	0.2775	0.0000	0.3613	0.8212	1.1908	2.3533	0.4801	2.1284	0.0999	0.4597	0.4801	0.4631	6.1526	0.2707
15006	ZCTA5 15006	MULTIPOLYGON (((-79.88531 4…	292	131	131	0	0	0	37	0	58	13	0	0	0	0	0	0	0	0	0	0.0	0.0	0.0	13.3	0.0	19.9	4.5	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0000	0.0000	0.0000	0.7911	0.0000	0.5091	0.0545	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000	0.7911	0.5636	0.0000	0.0000	0.0392	0.0148	0.0000	0.0000	1.3547	0.0136
15007	ZCTA5 15007	MULTIPOLYGON (((-79.93657 4…	629	167	167	387	0	51	0	19	414	22	295	0	0	0	0	0	0	0	0	61.5	0.0	30.5	0.0	3.0	65.8	3.5	46.9	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.9920	0.0000	0.8627	0.0000	0.3186	0.9881	0.0499	0.9903	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000	2.1733	2.0283	0.0000	0.0000	0.3785	0.2985	0.0000	0.0000	4.2016	0.1022
15009	ZCTA5 15009	MULTIPOLYGON (((-80.45619 4…	15114	7057	6509	1904	323	1424	404	315	3718	3253	2196	271	57	869	528	18	23	335	304	12.8	4.5	21.9	3.6	2.1	24.6	21.5	14.8	4.2	0.4	5.7	7.5	0.3	0.4	5.1	2.0	0.3109	0.5259	0.5436	0.1635	0.2186	0.7832	0.6317	0.4991	0.5322	0.6226	0.5176	0.8199	0.2479	0.4479	0.5105	0.8031	1.7625	3.0688	0.5176	2.8293	0.2236	0.8593	0.5176	0.6850	8.1782	0.6118
15010	ZCTA5 15010	MULTIPOLYGON (((-80.48268 4…	27205	12699	11363	4761	811	2322	1237	923	5815	5636	4584	820	42	3819	811	251	135	1243	1268	18.3	6.1	20.4	6.6	3.4	21.4	20.7	17.1	7.2	0.2	14.0	6.4	2.0	1.2	10.9	4.7	0.5359	0.7077	0.4536	0.3678	0.3759	0.6101	0.5647	0.6547	0.8114	0.5312	0.7611	0.7954	0.3806	0.6632	0.8046	0.9030	2.4409	3.1721	0.7611	3.5468	0.4784	0.8939	0.7611	0.9222	9.9209	0.8610
15012	ZCTA5 15012	MULTIPOLYGON (((-79.90417 4…	15243	7709	6883	2639	376	1186	638	470	3595	2750	2629	348	86	1048	353	445	22	481	144	17.3	5.0	17.2	5.7	3.1	23.6	18.0	17.2	5.1	0.6	6.9	4.6	5.8	0.3	7.0	0.9	0.4972	0.5806	0.2917	0.2991	0.3311	0.7412	0.3570	0.6581	0.6279	0.6924	0.5681	0.7356	0.5436	0.4137	0.6370	0.7003	1.9997	3.0766	0.5681	3.0302	0.3099	0.8627	0.5681	0.7491	8.6746	0.7015
15014	ZCTA5 15014	MULTIPOLYGON (((-79.74922 4…	3055	1671	1450	761	41	334	136	175	464	527	379	162	0	264	48	0	7	136	91	26.0	2.5	23.0	6.2	5.9	15.2	17.3	12.8	11.2	0.0	8.6	2.9	0.0	0.5	9.4	3.0	0.7625	0.2593	0.6057	0.3331	0.6599	0.2287	0.3121	0.3481	0.9242	0.0000	0.6328	0.6575	0.0000	0.4752	0.7533	0.8587	2.6205	1.8131	0.6328	2.7447	0.5573	0.1947	0.6328	0.6578	7.8111	0.5409

For visualization purposes, we’ll simplify the table and keep only the GEOID, geometry and SVIs.

svi <- pa_zcta_2020_geo_svi %>% 
  select(GEOID, geometry, contains("RPL_theme"))
svi %>% head(10) %>% kable()

GEOID	geometry	RPL_theme1	RPL_theme2	RPL_theme3	RPL_theme4	RPL_themes
15001	MULTIPOLYGON (((-80.43758 4…	0.5028	0.7418	0.8173	0.8224	0.8059
15003	MULTIPOLYGON (((-80.2366 40…	0.5380	0.8473	0.7860	0.6515	0.7707
15004	MULTIPOLYGON (((-80.38652 4…	0.4745	0.2520	0.6118	0.1266	0.2287
15005	MULTIPOLYGON (((-80.24356 4…	0.0999	0.4597	0.4801	0.4631	0.2707
15006	MULTIPOLYGON (((-79.88531 4…	0.0392	0.0148	0.0000	0.0000	0.0136
15007	MULTIPOLYGON (((-79.93657 4…	0.3785	0.2985	0.0000	0.0000	0.1022
15009	MULTIPOLYGON (((-80.45619 4…	0.2236	0.8593	0.5176	0.6850	0.6118
15010	MULTIPOLYGON (((-80.48268 4…	0.4784	0.8939	0.7611	0.9222	0.8610
15012	MULTIPOLYGON (((-79.90417 4…	0.3099	0.8627	0.5681	0.7491	0.7015
15014	MULTIPOLYGON (((-79.74922 4…	0.5573	0.1947	0.6328	0.6578	0.5409

Interactive maps for overall SVI

With the simple feature geometry, we can visualize SVI patterns and perform spatial SVI analysis with any mapping tool. Here, we’re using a powerful package leaflet for interactive maps.

First we’ll examine the missing value.

missing <- svi %>% filter(is.na(RPL_theme1)) 
missing %>% kable()

GEOID	geometry	RPL_theme1	RPL_theme2	RPL_theme3	RPL_theme4	RPL_themes
15260	MULTIPOLYGON (((-79.95539 4…	NA	NA	NA	NA	NA
15616	MULTIPOLYGON (((-79.563 40….	NA	NA	NA	NA	NA
15691	MULTIPOLYGON (((-79.68701 4…	NA	NA	NA	NA	NA
16312	MULTIPOLYGON (((-79.30513 4…	NA	NA	NA	NA	NA
17077	MULTIPOLYGON (((-76.53944 4…	NA	NA	NA	NA	NA
17120	MULTIPOLYGON (((-76.88651 4…	NA	NA	NA	NA	NA
17822	MULTIPOLYGON (((-76.60949 4…	NA	NA	NA	NA	NA
18225	MULTIPOLYGON (((-75.97559 4…	NA	NA	NA	NA	NA
18430	MULTIPOLYGON (((-75.49468 4…	NA	NA	NA	NA	NA
18936	MULTIPOLYGON (((-75.23686 4…	NA	NA	NA	NA	NA
19109	MULTIPOLYGON (((-75.16434 3…	NA	NA	NA	NA	NA
19112	MULTIPOLYGON (((-75.19692 3…	NA	NA	NA	NA	NA
19442	MULTIPOLYGON (((-75.58142 4…	NA	NA	NA	NA	NA

Looks like there are a few ZCTAs with missing values in SVIs. (Upon checking the full SVI table, we can see that most of the individual SVI variables are 0 for those ZCTAs.) So we’ll remove those ZCTAs for better visualization.

To set up the interactive map for overall SVI, leaflet() does most of the heavy lifting. Here, we’ll just add some customized color palette and labels for the appearance.

svi_clean <- svi %>% drop_na()

#set color palette
pal <- colorNumeric(
  palette = c("orange","navy"),
  domain = svi_clean$RPL_themes
)

#set label
zcta_label <- glue("<h3 style='margin: 0px'>{svi_clean$GEOID}</h3>
                    overall SVI: {svi_clean$RPL_themes}") %>%
  map(~HTML(.x))

leaflet(svi_clean) %>% 
  addProviderTiles(providers$CartoDB.Voyager) %>% 
  addPolygons(color = "white", 
              weight = 0.5,
              smoothFactor = 0.5,
              opacity = 1,
              fillColor = ~pal(RPL_themes),
              fillOpacity = 0.8,
              highlightOptions = highlightOptions(
                weight = 5,
                color = "white",
                fillOpacity = 0.8,
                bringToFront = TRUE),
              label = zcta_label,
              labelOptions = labelOptions(
                style = list(
                  "font-family" = "Fira Sans, sans-serif",
                  "font-size" = "1.2em"
                ))
              ) %>% 
  addLegend("bottomleft",
            pal = pal,
            values = ~RPL_themes,
            title = "Overall SVI in all ZCTAs in PA (2020)",
            #labFormat = labelFormat(prefix = "$"),
            opacity = 1)
#> Warning: sf layer has inconsistent datum (+proj=longlat +datum=NAD83 +no_defs).
#> Need '+proj=longlat +datum=WGS84'

From the interactive map, we can visualize easily how SVI varies in different regions in PA and zoom in to examine specific ZCTAs of interest, making it a helpful approach to explore new ideas, patterns and analyses.